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ABSTRACT 
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We show that condensation is an efficient particle growth mechanism, leading to growth beyond decimeter-sized pebbles close to 
an ice line in protoplanetary discs. As coagulation of dust particles is frustrated by bouncing and fragmentation, condensation could 
be a complementary, or even dominant, growth mode in the early stages of planet formation. Ice particles diffuse across the ice line 
and sublimate, and vapour diffusing back across the ice line recondenses onto already existing particles, causing them to grow. We 
develop a numerical model of the dynamical behaviour of ice particles close to the water ice line, approximately 3 AU from the host 
star. Particles move with the turbulent gas, modelled as a random walk. They also sediment towards the midplane and drift radially 
towards the central star. Condensation and sublimation are calculated using a Monte Carlo approach. Our results indicate that, with 
a turbulent o--value of 0.01, growth from millimeter to at least decimeter-sized pebbles is possible on a time scale of 1000 years. We 
find that particle growth is dominated by ice and vapour transport across the radial ice line, with growth due to transport across the 
atmospheric ice line being negligible. Ice particles mix outwards by turbulent diffusion, leading to net growth across the entire cold 
region. The resulting particles are large enough to be sensitive to concentration by streaming instabilities, and in pressure bumps and 
vortices, which can cause further growth into planetesimals. In our model, particles are considered to be homogeneous ice particles. 
Taking into account the more realistic composition of ice condensed onto rocky ice nuclei might affect the growth time scales, by 
release of refractory ice nuclei after sublimation. We also ignore sticking and fragmentation in particle collisions. These effects will 
be the subject of future investigations. 
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1. Introduction 

Planets form in protoplanetary discs of gas and dust, surround- 
ing young stars. In the classical planet formation scenario dust 
grains coll i de, st ick together and form larger and larger bodies 
(Saf ronovl[T969h . Particles of up to millimeter-sizes stick to- 
gether due to contact forces and kilometer-sized and larger bod- 
ies are held together by gravity. How particles grow from mil- 
limeter to kilometer-sized planetesimals is difficult to explain 
by collisions, as these particles tend to fragme nt or bounce off 
each other as they collide, instead of sticking (Blu m & Wurml 
2008). Further, solids on Keplerian orbits experience a head- 
wind from the slow-orbiting pressure-supported gas, and there- 
fore lose angular momentu m and drift in towards the central 
star dWeidenschillingl Il977l) . This radial drift velocity peaks for 
meter-sized particles, which drift into the star fr om 1 AU on a 
time-scale of only 100 orbits dBrauer et al.ll2008l) . 

Once particles grow to centimeter-sized to decimeter-sized 
pebbles, further growth into planetesimals is possible via parti- 
cle concentration and gravitational collapse. The streaming in- 
stability causes particles to clump together due to the velocity 
difference between the particles and the gas in the disc, and 
planetesimals form t h rough subsequent grav i tation al collapse 
(lYoudin & Goodman! 120051: I Johansen et ail 120071). Particles 
also c oncentrate in long-live d vorti ces dBarge & S ommeria, 
[1991 iKlahr & Bodenheimel 12003b and pressure bumps 
dJohansen et all l2009t) excited in the turbulent flow. However, 
an efficient mechanism for growth from millimeter-sized dust 
grains to pebbles is needed to explain the formation of particles 
that are large enough to take part in such concentration events. 



Sticking, or coagulation, as a growth mechanism has 



been extensively s tudied, both experimentally (Blum & Wurm, 
2008 :IZsom et all 120101) and numerically dBraueret all 12008: 
Windmark et al However, an often overlooked concept 

in planet form ation models is that of gro wth via condensation 
near ice lines. IStevenson & Lunina (1 19881) investigated material 
enhancement at the ice line due to diffusive redistribution and 
subsequent condensation of water vapour from the inner part of 
the disc. With the assumptions of efficient condensation by ho- 
mogeneous nucleation, and by ignoring inwards transport of ma- 
terial, the a uthors found a significa nt solid density enhancement. 
In contrast, ICuzzi & Zahnld (120041) considered inwards material 
drift across the ice line, which would enhance the vapour density 
in the inner disc, followed by a phase of accumulation of solids 
by condensation onto immobile planetesimals formed outside of 
the ice line. 

In this work the dynamical behavior and growth of small 
(~1 mm) seed particles via ice condensation is studied in com- 
puter simulations. Ice particles moving with the turbulent gas 
encounter the water ice line where ice sublimates and becomes 
water vapour. We consider both the radial ice line, separating 
the inner, hotter part of the disc from the colder region further 
away from the star, and the atmospheric ice line, which sepa- 
rates the hot atmosphere of the disc from the cold midp lane layer 
dDuUemond & DominikL 12004b iMeiierink etafl 120091) . Due to 
the turbulent motion of the gas in the protoplanetary disc wa- 
ter vapour diffuses across the ice line into the condensation re- 
gion where it condenses onto existing ice particles. Many ice 
particles will move across the ice line and sublimate. However, 
since small particles are coupled to the gas and thereby effec- 
tively move in a random walk due to turbulence, a significant 
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fraction will be lucky enough to stay within the condensation 
region of the disc, growing to larger sizes as diffusing water 
vapour condenses onto them. The total mass in the ice compo- 
nent is maintained during this process, as the loss of ice particles 
to sublimation is completely compensated for by the growth of 
the remaining particles. 

In a moderately turbulent disc we find particle growth to be- 
yond decimeter-sizes on a time-scale of a few 1000 orbits. The 
most significant growth takes place locally, close to the radial 
ice line. However, efficient radial mixing supplies the midplane 
with large particles even several scale heights outside of the ice 
line. The atmospheric ice line, located at z > 3 H where disc ma- 
terial is heated directly by the central star (Chiang & Gol dreichL 
Il997l) . contributes less to the particle growth, generally on the 
order of jjm. 

The radial water ice line, or c ondensation front, is found at 
r ~ 3 AU around a solar type star dLecar et al. , 2006). However, 
growth by condensation is not only applicable to water ice, but 
also to any other volatile found in the protoplanetary disc. Other 
important condensation fronts relevant for planet formation in- 
clude those of ammonia, methane, carbon monoxide and molec- 
ular nitrogen at much la rger radii than the water ice li ne, and 
of silicates at r « 1 AU dLoddersl 12001 lOi et all [201 lb . These 
condensation fronts will be the topic of a future study. 

The paper is organised as follows. In Section 2 the numerical 
model used for particle dynamics and growth by condensation 
is described, including the units used in this paper. The Monte 
Carlo scheme for condensation is further explained in Section 
3. In Section 4 we describe the test problems used to validate 
the code. We present our results in Section 5, and investigate the 
influence of the position of the atmospheric ice line, the turbu- 
lence strength and the effects of radial mixing. In a moderately 
turbulent disc (a = 10~ 2 ), pebbles of centimeters to decimeters 
form within a few 1000 years. We discuss general assumptions 
and simplifications made in Section 6. Finally, in Section 7 we 
conclude that particle growth by condensation is an important 
mechanism for forming pebbles in protoplanetary discs, which 
should be taken into account in future studies of planetesimal 
formation. 



2. Numerical model 

2.1. Simulation box, units and boundary conditions 

We simulate a two-dimensional region around the water ice line. 
The simulation domain is set in the radial r and vertical z di- 
rection, whereas the azimuthal direction is ignored, for simplic- 
ity and under the assumption of axisymmetry. The fundamental 
units are the gas scale height H as a length unit, sound speed 
c s as a velocity unit and inverse Keplerian orbits Q l as a time 
unit. The relation between these quantities is H = c s /Q. Setting 
H — c s = Q — 1 gives a system which is scalable to any region 
of the protoplanetary disc. Hence our results apply equally well 
to other condensation fronts, such as those of CO or silicates, al- 
though the transition and scaling between drag regimes and the 
fractional ice abundance depend on the choice of radial location 
in the disc. Here we interpret the results in terms of the water ice 
line at arou nd 3 AU from the star, using a water ice mass fraction 
of 0.57 1 % dLoddersl 120031) . 

Condensation is considered as a neighbour interaction. A lin- 
ear scaling between number of particles and number of calcula- 
tions is obtained by use of a mapping scheme, where the simula- 
tion domain is divided into a number of grid cells, with particle 



interactions possible between particles within the same grid cell 
only. 

Initially, both vapour and ice particles have a Gaussian dis- 
tribution in the vertical direction, following the gas distribution. 
At the beginning of a simulation particles are set to be ice or 
vapour, depending on whether they are located within or outside 
of the condensation region. 

2.2. Particle sizes 

Particles are coupled to the gas via drag forces and can be char- 
acteri zed by their dimensionless friction time (IWeidensch illing. 
fl977h . or Stokes number, 



(Ep) 



(St) 



Rp. 

4 R 2 p. 

9Hp g A 



for R < (9/4) A, 



for R > (9/ 4) A, 



(1) 



(2) 



where A is the mean free path of the gas and p. and p g are the 
material and gas density. The superscripts (Ep) and (St) denote 
Epstein and Stokes drag regime, the two drag regimes relevant 
for the particle sizes and gas density considered in this paper. 
We formulate the particle radius R in units of Rj , where R/Ri — 



Qrf p) so that 



#1 



P. 



(3) 



To find the corresponding particle size in meters, we use the ma- 
terial density of ice p. as 1 g cirT 3 , the midplane gas density p„ 
and column densit y £ from the M inimum Mass Solar Nebula 
(MMSN) model of lHavashil fl98lh . 



Pg = 



2T(r) 



y/2nH(r) ' 



^noogcnr^)- 



1.5 



(4) 
(5) 

(6) 



giving the scaling of 

Ri ~ 1.3mf- — I 

\3AU/ 

close to the water ice line. 
2.3. Condensation scheme 



The rate of change in particle mass due to condensation and sub- 
limation is 



dm . n2 L Psat 



(7) 



dSupulver &LinL H000). Here, Vth is the thermal velocity of 
vapour, p v is the vapour density and P sat and P v is the satu- 
rated vapour pressure and vapour pressure, respectively. Both 
the vapour pressure and the saturated vapour pressure can be ex- 
pressed as densities using the ideal gas law. The vapour pressure 
can be written as 



k B T 
P v = Pv • 

m v 



(8) 



and the saturated vapour pressure can be expressed in an equiv- 
alent way. Here, k% is the Boltzmann constant, T is the temper- 
ature, p v is the vapour density and m v is the mass of one vapour 
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particle. Assuming spherical particles Eq. |7]can be rewritten in 
terms of particle radius R and using vapour densities instead of 
pressures, 



dR vth , . 

~T~ — (Pv - Psat) • 

at p. 



(9) 



In the limit where p v <K p sat sublimation dominates the time 
evolution. The sublimation time scale in terms of particle radius 
is then 



and in terms of mass 



Rp. 

VthPsat 



1 Rp. 



3 V t hp sa t 



(10) 
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Correspondingly, when p sat «: p v condensation dominates, with 
condensation time scale in terms of radius 



T C ,R = 



Rp. 

VthPv 



and in terms of mass 



1 Rp. 

Tr = -- 



3 v th p v 



(12) 



(13) 



In equilibrium p v = p sat with t s = t c . By definition p v = fn^oP°, 
where fu 2 o is the water mass fraction, at an ice line. The rapid 
increase of the temperature into the condensation region implies 
Psat » Pg and hence t s <k Tf , allowing for modelling of sublima- 
tion as an instantaneous process as long as the particle friction 
time, Tf, is shorter than the characteristic time-scale of the tur- 
bulent eddies, Or 1 . 

Condensation can not be assumed to occur instantaneously, 
as p v never increases much beyond /h,oPg- The condensation 
time scale is proportional to the size of the particle, with the 
dimensionless condensation time scale being 



Or 



(Ep) 



(14) 



Vth, where pn,o and p^ are the mean 



assuming c s « 

molecular weights of water and molecular hydrogen. The con- 
densation process is modelled by a Monte Carlo approach, where 
the probability of condensation for vapour onto an ice particle in 
the condensation zone is proportional to the size of the available 
ice surface. This is further described in Section [3] 

2.4. Superparticle approach 

The ice and vapour components are modelled using a su- 
perparticle approach , relat ed to the coagulation algorithm of 
Zsom & Dullemond (2008), where a superparticle is a numer- 
ical representation of a large number of physical particles with 
identical properties. Here the most important properties are the 
physical state (ice or vapour), size of constituent particles if ice, 
and material density. The number of superparticles N is much 
smaller than the number of physical particles, so that superpar- 
ticle i represents n, physical particles. In these simulations, typ- 
ically N = 10000 for a simulation area of 12 H in the vertical 
direction and 1 .5 H in the radial direction, a number chosen to be 
computationally easy to handle. The number of superparticles is 
fixed throughout a simulation, and each superparticle represents 




Fig. 1. Sketch of sublimation and condensation for superparti- 
cles. Blue represents ice and red vapour. Left panel: Ice super- 
particle, representing a mass M in the form of 4 physical ice 
particles, each of mass m, sublimates and turns into a vapour su- 
perparticle of mass M. Right panel: Vapour superparticle of mass 
M condenses onto an ice superparticle of mass M, representing 
the mass 4 m. The physical ice particles are after the event repre- 
sented by two ice superparticles, each representing a total mass 
M, but now in the form of 2 physical ice particles each with mass 
2 m. Both mass and the number of superparticles are conserved 
between fo and t\ . In the case of condensation, also the number 
of physical particles is conserved. 



an equal and fixed mass M that does not change. The mass m 
of the physical particles represented by the superparticle does 
however change, and so does n,-, the number of physical parti- 
cles represented by a superparticle, as the total mass represented 
by the superparticle is always M = ;«,■ n,. Also when an ice par- 
ticle undergoes a phase transition and becomes vapour the total 
mass M of the superparticle stays constant. The superparticle 
approach to condensation and sublimation is sketched in Fig.Q] 

2.5. Random walk of particles 

In this section the modelling of particle motions is described. 
We start by introducing the random walk used to describe the 
turbulent motions of the gas, sufficient for modelling small par- 
ticle motions, and then move on to include the more complex 
behaviour of larger particles. As particles grow larger, additional 
effects need to be taken into account. Firstly, the random step 
length becomes smaller, as larger particles follow the turbulent 
eddies a shorter length each time step. Secondly, for larger par- 
ticles gravity towards the midplane has an increasingly impor- 
tant effect, causing particles to sediment towards the midplane. 
Finally, radial drift towards the central star must be taken into 
account. 

With the dim e nsionl ess constant a introduced by 
Shakura & Sunvaev] (Il973l) . the turbulent diffusion coeffi- 
cient of gas and small particles in an accretion disc can be 
written as 

D = ac s H. (15) 

We here assume that particles move with velocity v ~ -\Jac s , and 
that the length a particle moves coherently is / ~ -\/aH, which 
gives rise to D ~ vl ~ ac s H by mixing length arguments. The 
correlation time is defined as 



v 



(16) 



This means that the effect of turbulence on a small parti- 
cle is to move it in a random direction with a characteristic 
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speed v, during the correlation time r. The characteristic length- 
scale, velocity-scale and time-scale are all set by the turbu- 
lence. We use a default a-value of 10~ 2 , motivated by observa - 
tions of accretion rates in young stars (Hartma nn et al.l [l998). 
Howe ver, the prescence of a dead zone can give a low er a- 
value dFleming & StoneL 120031: lOishi & Mac Lowi [2009), and 
there are al so observations suggesting a = 10~ 4 in protoplan- 
etary discs dMulders & Dorruifikl |2012|) . Therefore we also run 
simulations with a = 10~ 3 and a = 10~ 4 . 

The distance a particle moves each time step, where one time 
step is At = Q l , is found by equating the generic expression for 
the diffusion coefficient of a random walk in 2D, D = I 2 1 (At), 
with that of diffusion described by Eq. Q3] giving 



I = 2 Vatf • 



(17) 



Practically, the random walk for particles coupled to the gas is 
modelled by setting the length the particle moves in the radial 
direction Ar and in the vertical direction Az to 



(Ar)™, = I cos 9 , 
(Az)™ = I sin 6 . 



(18) 
(19) 



Here < 8 < 2n is chosen randomly for each particle and time 
step, giving a random direction, but a fixed length, for how a par- 
ticle moves in one time step. This forms the basis of the particle 
dynamics used in this code, and is valid for small particles per- 
fectly coupled to the turbulent gas. However, for larger particles 
this simple approximation breaks down. A number of corrections 
to the random walk model will now be introduced, one at a time, 
before arriving at the final model for the particle dynamics. 

Firstly, the random step length must be adjusted, as larger 
particles move shorter distances with the turbulent eddies each 
time step. The diffusion coefficient D = vl = I 2 It = / 2 Q is 
modified to 

° (20) 



D 



1 + (Qt s ) 2 



following the ana ly tical theory developed and tested by 
lYoudin & Lithwickl d2007h . This gives a modified size- 
dependent step length of 



r 



/ 



2y/aH 

Vl + (QT f ) 2 Vl + ( Q ff) 2 



(21) 



For small particles I* as I, whereas for larger particles P < I. 
Separating the radial and vertical directions, Eqs. ITS"! and [T91 are 
replaced by 



(Ar) r 
(Az) r 



/ cos 6 



Vl + (^f) 2 
/ sin 

Vi + (^f) 2 



(22) 
(23) 



As particles grow larger, and drag forces decrease in strength, 
the gravity towards the midplane influences the particle motion 
more, so that large particles sediment to the midplane. Since this 
affects the vertical direction only, it is only necessary to modify 
Az. This is done by changing the step size in the vertical direc- 
tion to (Az) tot = (Az)™ + (Az) sed , where (Az) sed is an increasing 
function of particle size. The sedimentation step length is set by 
the terminal velocity of the particles \\ = -TfQ 2 z. However, in 
order to modify the step length correctly, we must ensure that 
the resulting particle scale height is the same as the analytically 
expected scale height. The analytically expected particle scale 



height, for a given particle size and turbulence strength, result- 
ing from diffusion and sedimentation is 



H 



Qrt + a 



(24) 



as determined in compu ter simulations of part icles in turbu- 
lent protoplanetary discs (Carballi do et al.l [2006b and explained 
in the analytical diffusion framework of lYoudin & Lithwickl 
J2001 . This can now be compared with the particle scale height 
arising from considering a combination of a random walk and 
terminal velocity sedimentation, which is what is used for mod- 
elling. The random walk of particles mimics turbulent diffusion 
with coefficient D = ac s H, and sedimentation-diffusion equilib- 
rium occurs at particle scale height 



Hi 
H 



Qr~ f 



(25) 



To construct a model that gives the correct scale height, as given 
by Eq. [24] we therefore need to modify the dimensionless fric- 
tion time in the terminal velocity expression to 



(QTf)* = Qt{ + a . 



(26) 



Effectively we let even tiny particles sediment in order to fol- 
low the Gaussian distribution of the gas, rather than the uni- 
form distribution resulting from a pure random walk. Changing 
QTf — > (QTf)* in Eq. [25] reproduces the correct particle scale 
height in Eq. [24] With the modified dimensionless friction time 
(QTf)* in the terminal velocity expression, the step size in the 
vertical direction can be written as 



(Az) to t = (Az)™ - (QTf + a)QzAt . 



(27) 



However, this expression leads to a numerical instability in 
which large particles overshoot the midplane as they sediment, 
since for QTf > 1 and At = Q~ l we get 



(QTffQzAt > z . 



(28) 



This gives particle oscillations that are amplified in time. To 
eliminate this problem, Eq. [27]is modified to 



(Az) tot = (Az)nv - 



(QTf + a)QzAt 

1 + (QTf) 2 



(29) 



following the scheme o f Youdi n & Lithwickl d2007l) . which elim- 
inates the amplifying of the vertical particle oscillations. Instead 
it gives larger particles a longer settling time, corresponding to 
the time they would have spent oscillating before settling. For a 
comparison between the particle scale height resulting from Eq. 
[29] and the theoretical scale height from Eq. [24] see Fig. [2] The 
correlation between the analytical and modelled curve in this fig- 
ure shows that the particle dynamics in the vertical direction are 
correctly modelled by Eq. [29] 

The step in the radial direction needs to be modified due to 
the radial drift towards the central star. This is also particle size 
dependent, and (Ar) tot is modified to 



(Ar) tot = (Ar)™ + (Ar) rd = (Ar) r 



2??vk 



QTf + (QTfY 



■At, (30) 



where Vr is the Keplerian velocity and 77VK is the veloc- 
ity difference b e tween the gas and the particles, following 
Weidenschillin 2 (fl977h . This gives a radial drift velocity that 
peaks for QTf = 1 particles. We use 77VK = 0.05c s , a 
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reasonable value at 3A U in the MMSN (ICuzzi et all Il993t 
IChiang& Youdinll2010h . 

The friction time of a particle is inversely proportional to the 
gas density, QTf oc p" 1 , which decreases in the vertical direc- 
tion with distance from the midplane. Taking into account the 
Gaussian distribution of the gas around the midplane, the dimen- 
sionless friction time of a particle is modified according to 



Qrf — > QTf e , 



where 



exp 



[z 2 /(2H 2 )] 



(31) 



(32) 



Simulations have shown that the heterogeneous vertical gas dis- 
tribution leads to a turbulent a- value increa sing with height 
above the midplane (Fromang & Nelson! 120091) . We do however 
not take this effect into account in our simulations. 

For large simulation domains radial gas density gradients are 
also of importance, and the random walk would then have to be 
adjusted in order to mimic th e motion of especially smal l parti- 
cles, well coupled to the gas (Hug hes & Armitagell2010h . Since 
the model used in this paper is local, we ignore this effect. 

The equations used to describe particle dynamics in the code 
are thus 



Ar = 



/ sin6» 



(QT f e + a)Qz 
Vl + (QTf£) 2 1 + ( Qt ( £ ) 2 



At, 



Vl + (QT t e) 2 eQTf + (eQTfY 



-At , 



(33) 
(34) 



where the first term for both the radial and vertical direction de- 
scribe the random walk, with a step length modified according 
to particle size, and the second term describes sedimentation and 
radial drift, for the vertical and radial direction, respectively. 

2.6. Transition between drag regimes 

For small particles with a radius of less than (9/4) of the gas 
mean free path, Epstein drag applies, whereas for larger particle 
Stokes drag is the relevant regime. The transition size given by 
Eqs. [T]and [2] can be rewritten in dimensionless friction time as 



(QTf) 



(Ep) 



9 p. A 9tt p.fxH 



4p s H 2 2?o- a 



0.5. 



(35) 



when applied to the midplane of the disc. The mean molecu- 
lar weight is p — 3.9 x 10~ 24 g, the molecula r cross-section of 
molec ular hydrogen is cr mo \ — 2 x lCT 15 cm 2 (Nakagawa et all 
1986) and the column density Z at 3 AU is given by Eq.|5] 

Converting dimensionless friction time in the Epstein drag 
regime to the Stokes drag regime is done using 



Qt\ 



.(St) 



2 E cr mol 
9n p.pH 



(Q T f p) f * 2.0 (Q T f p) f , (36) 



where the scale height is given by the MMSN model (Havashi, 
119811) as 



H = 0.033 



1.25 



(37) 



3. Condensation in a Monte Carlo scheme 

We treat condensation as a neighbour interaction, and require 
vapour to be present near an ice particle for interaction to oc- 
cur. In this model "near" means being in the same grid cell, so 
that it is possible for an ice particle to interact only with vapour 
within the same grid cell. Condensation in each time step is thus 
decided one grid cell at a time. Any grid cell can, and most of- 
ten does, harbour more than one ice particle. Therefore, whether 
or not condensation takes place in a given grid cell and a given 
time step, is decided in a two-step process for each vapour parti- 
cle present in the grid cell. 

The two-step procedure is here described for one grid cell 
and one time step At, with one vapour particle present in this 
grid cell. In the case where more than one vapour particle is 
present in the grid cell, the two-step procedure is repeated for 
each vapour particle present. The first step is to decide whether 
this vapour particle condenses onto any of the ice particles. If 
it does, the second step is to decide which of the ice particles 
it condenses onto. In this scheme, any vapour particle can only 
condense onto one single ice particle, so it can not be shared 
amongst several particles. One ice particle can on the other hand 
experience several growth events during one time step, if several 
vapour particles are present. 

As a first step it is decided whether or not the vapour particle 
is involved in a condensation event during the time step At. In 
order to do this the total interaction probability for all ice parti- 
cles in the grid cell is needed. The interaction probability for ice 
particle i is 

Pi = l-exp(-Ar/T;), (38) 

where t,- is the condensation time scale for particle i, given by 
Eq. [14] The expression for interaction probability is chosen to 
always give </>,■< 1, for any Af, also for At » r,. The algo- 
rithm proceeds by calculating the probability that no interaction 
occurs. For n ice particles, the total probability that no interac- 
tion occurs is 



Po = (1 - Pi) ■ (1 - pi) ■ 



d-Pn). 



(39) 



To decide whether condensation occurs, a random number r\ e 
[0, 1 [ is generated. If r\ is smaller than po nothing happens, 
whereas if r\ is larger than po vapour condenses onto one of 
the ice particles. 

If condensation occurs, a second step is needed in order to 
decide onto which of the ice particles in the grid cell. In this step 
it is no longer the absolute probability that is of interest, but the 
relative probabilities of the ice particles in the grid cell. This can 
be expressed as 

Af 

P* = - ■ (40) 

The relative probability for each of the ice particles in the grid 
cell are placed in a sequence, 



0, p\, p\ +p* 2 , p\ + p* 2 +p* 3 , 



(41) 



and a new random number is generated, r 2 6 [0, 2; P*A ■ Which 
interval in the sequence r 2 falls in decides which ice particle the 
vapour particle condenses onto, such that r 2 falling in the interval 
]p*_ l ,p*] implies condensation onto ice particle i. 

The condensation time scale in Eq.[13]is proportional to the 
radius of the particle, so that smaller particles have a shorter, 
and larger particles a longer, condensation time scale. This cor- 
responds to a larger condensation probability for smaller par- 
ticles and vice versa, which may seem counter-intuitive, but 
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is explained by the fact that the simulation handles superpar- 
ticles, each representing the same mass M. Therefore a super- 
particle representing small physical particles represents a larger 
combined surface area than one representing large particles. 
Condensation is thus more likely to happen to a superparti- 
cle representing small particles. Looking at a single-physical- 
particle-basis, the result is correct averaged over many timesteps 
and particles. A condensation event means doubling the mass of 
the physical particles involved, implying that a small particle ex- 
periences many condensation events, but the absolute growth in 
radius each event is small, whereas a single large physical parti- 
cle experiences few condensation events, but with a large growth 
in absolute radius each event. 

4. Test problems 

In order to understand the results of the computer simulations 
and to make sure that the code is functioning correctly, two ma- 
jor tests of the algorithms used were made. Firstly, the dynam- 
ical behaviour of the code was tested without including parti- 
cle growth. Secondly, the particle growth algorithm was tested, 
without taking spatial dimensions into account. 

4.1. Test of the dynamical behaviour of the particles 

The dynamical behaviour of the particles is tested by excluding 
condensation and sublimation from the simulation. This means 
that particles move due to turbulence, stirring them up in a turbu- 
lent diffusion, and due to gravity directed towards the midplane, 
but no particle growth is included. We also ignore radial drift. 
The particles settle to an equilibrium, where the particle scale 
height depends on the size of the particles, since the strength of 
the coupling to the turbulent gas is a function of particle size. 
Fig. [2] shows the particle scale height relative to the gas scale 
height Hp/H as a function of particle size R. The particles settle 
into a Gaussian distribution around the midplane, so the parti- 
cle scale height in equilibrium can be retrieved as the root mean 
square of the particle positions, 




where n is the number of particles and z, is the vertical position 
of particles. This is compared to the theoretically expected parti- 
cle scale height, given by an equilibrium between sedimentation 
and turbulent diffusion, 

Fig. [2] shows how the particle scale height depends on the 
size of the particles, with the modelled values represented by the 
black full line and the expected analytical values as a red dot- 
ted line. The modelled line clearly follows the analytical curve, 
showing that the dynamical behaviour of the particles is as ex- 
pected. The particle size is given in units of Ri, where fiixlm 
at r — 3 AU. From the figure one can see that small parti- 
cles are fully coupled to the gas, as H v /H * 1 for particles of 
R < 10 _3 i?i, corresponding to mm-sized ice particles near the 
water ice line. Further, a slight change of the slope can be seen 
at R ss 0.5 R\. This is caused by the change from Epstein to 
Stokes drag regime, which makes the particles decouple more 
quickly from the gas and hence increases sedimentation. 




Fig. 2. Comparison between modelled and analytical particle 
scale height as a function of particle size, for a turbulence 
strength of a — 10~ 2 . The black full line denotes modelled 
values and the red dotted line analytical values. Small particles 
R < 1CT 3 are perfectly coupled to the gas, whereas larger parti- 
cles sediment towards the midplane. The slope change at R » 0.5 
is due to the change from Epstein to Stokes drag regime. The par- 
ticle size is given in units of R\, where R\ * lmatr = 3 AU. 
Overall, modelled and analytical values are in good agreement. 



4.2. Test of the condensation algorithm 

We test the condensation algorithm by letting particles grow via 
condensation, without including the spatial dimensions. The ice 
line is therefore not included, so sublimation is not taken into 
account. The number of ice superparticles is known, as is the 
number of vapour superparticles. The total mass available for 
growth is known, since we run the simulation until all vapour 
has condensed onto the ice particles. By tracking the initial and 
final mass of all superparticles the growth in radius is followed, 
and can be compared to a theoretical expectation. 

In reality, condensation is a continuous process on macro- 
scopic scales, where single water molecules condense onto ice 
particles so that growth happens little by little all the time. In the 
Monte Carlo scheme used in this code this continuous behaviour 
is modelled by a discrete growth process, where a condensation 
event involves one ice superparticle and one vapour superpar- 
ticle. The mass of the physical ice particles represented by the 
ice superparticle doubles as the vapour superparticle completely 
condenses onto the ice, implying that there is no vapour super- 
particle left after the event. To conserve the number of ice parti- 
cles and superparticles (for easier coding purposes) this implies 
that the physical ice particles that before the event were repre- 
sented by the one ice superparticle involved in the event, now are 
represented by two superparticles, i.e. the number of physical ice 
particles before the event are after the event split between two 
ice superparticles. The mass of the physical particles involved in 
an event thus doubles as a condensation event occurs. As parti- 
cle radius is connected to mass as m oc R\ for any number of 
condensation events n ce the radius increases as 

R _> (2"«) 1/3 fl. (44) 
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Fig. 3. Comparison between average modelled growth and ana- 
lytically expected growth of ice particles. The final particle size 
as a function of initial particle size is shown for 2000 ice super- 
particles distributed evenly in 20 different initial size bins over 
the range R/R\ = [0, 1] , and with 20000 vapour superparticles 
added. Black crosses show modelled values and blue stars show 
the average growth in each size bin. The modelled average size 
roughly agrees with the analytical value, shown as a blue line. 
Red lines denote expected size for a number of condensation 
events between (lower line) and 10 (upper line). 

The final radius Rf\ m i is plotted as a function of R m i t as black 
crosses in Fig. [3] In this test 2000 ice superparticles were dis- 
tributed evenly in 20 different initial size bins over the range 
R/R\ = [0, 1] , and 20000 vapour superparticles were added. 
The possible growth in radius from mass doubling events is plot- 
ted as red lines in Fig. [3] for a different number of condensation 
events < n ce < 10. From the figure it is clear that all modelled 
values indeed fall on these lines, and never in-between, as ex- 
pected because of the discrete nature of the Monte Carlo method 
used. 

The modelled particle growth can be compared to the ana- 
lytically expected AR. For a physical ice particle i the total mass 
growth due to condensation during the time At is 



Am; 



(45) 



where p, is the material density of ice, Rqj is the initial and R, 
the final radius of the ice particle. However, the physical parti- 
cles are represented by superparticles in the code, and the mass 
growth must therefore be given for superparticles and not phys- 
ical particles. One superparticle represents the mass M, = «,/«,, 
where n, is the number of physical particles represented by the 
superparticles. 

The total mass of vapour that condenses onto the ice particles 
is given by the difference between the sum of final and initial ice 
superparticles, 



M v 



4tt 

: PvV = —p. 



(46) 



ice superparticles is not. As the physical particles represented by 
one superparticle grow in mass, the number of physical particles 
represented by the superparticle decreases since the total mass 
of a superparticle is constant. In order to conserve the number 
of physical particles, total number of superparticles and the to- 
tal mass of each superparticle, the number of ice superparticles 
increases as the physical particles grow in mass, so that each su- 
perparticle represents fewer and fewer physical particles , which 
is why the two sums in Eq.[46]are over different numbers of su- 
perparticles. 

The right handside of Eq. |46]can be rewritten using 



n = ^- = 

' nij 



PiV 
P.-Rj 



(47) 



Here V is the volume represented by a superparticle, which in 
this one-grid-cell test problem is equal to the total volume of the 
box. By cancelling terms and gathering the densities on the left 
handside Eq. [46]can be rewritten as 



(48) 



where N p and N p $ is the final and initial number of ice superpar- 
ticles, respectively. All superparticles represent the same density, 
and therefore the ratio between the two densities corresponds to 
the ratio between the number of vapour and ice superparticles 
represented, so that 

Pl = ^l, (49) 
Pi 1 

with N v being the number of vapour superparticles added to the 
simulation. Eg. l46l therefore is equivalent to the very simple ex- 
pression 

N v = N p - A/p, , (50) 

which the code can be checked against, as both the number of 
vapour superparticles put into the system, and the initial number 
of ice superparticles, is specified, and the final number of super- 
particles is given as an output from the program. This analysis 
shows that the code is mass-conserving, by construction. 

From Eq. [46] the theoretically expected AR, equal for all 
particles, can be found. The final radius can be rewritten as 
a, = AR + Rqj, where Roj is the initial radius of the particle. 
The two sums are not necessarily over the same number of su- 
perparticles, as each condensation event introduces a new ice 
superparticle in the system, at the expense of the vapour parti- 
cle representing the condensing vapour. The physical particles 
originally represented by the first ice superparticle are now split 
between the new and the old ice superparticle. The initial size of 
a physical particle represented by a converted vapour superpar- 
ticle can therefore be found as the initial size represented by the 
first ice superparticle. Eq. |46]is thereby rewritten as 



An 

M v = —p. 



J^miAR + Rorf-YjURli 

i 0,i 



(51) 



which, by using Eq. |47] cancelling terms and noting that (as in 
Eq.[46]lp v /p, = Nv, gives 



(52) 



The total number of physical ice particles is constant, as is the 
total number of ice and vapour superparticles, but the number of 



The initial number of vapour and ice superparticles N v and N p $ 
are known, as they are input parameters to the simulation. The 
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Table 1. List of simulations, with references to relevant figures in the second column. The dimensions of the simulation box are 
given as r m i n , r max , Zmm and z max in gas scale heights H. The number of grid cells in the radial and vertical direction is given as 
n r xn z and the number of particles as n p . Position of ice lines is given as Zi ce / H for the atmospheric ice line, and r; ce / H for the radial 
ice line when applicable. Turbulence strength is given by a and the initial particle size /?j n ; t is given in units of R\, approximately 
equivalent to size in meters. 



Run 


Fig. 


r,mn/H 




Zmm/H 


Zmax/H 


n, x n z 
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a 
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2a 
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2b 
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io- 3 


4a 


LlU 


-0.75 


0.75 


-6.0 


6.0 


64x512 


10 4 
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IO" 1 . 


io- 


-\ io- 1 


io- 3 


4b 


EES 


-0.75 


2.25 


-6.0 


6.0 


128 x 512 


2x 10 4 


3.0 


-0.3 


io- 4 , 


io- 


-\ io- 2 


IO" 3 


4c 


mi 


-0.75 


3.75 


-6.0 


6.0 


192 x 512 


3 x 10 4 


3.0 


-0.3 


10" 4 , 


io- 
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radius represented by a superparticle in the end of the simula- 
tion, ai, and when it is created, R c j, are both given as output 
from the code. The radius growth for each particle AR can there- 
fore be found. For N v = 20000 and 7Y p , = 2000, a value of 
AR m 0.202 is found. In Fig. [3] the blue full line represents the 
expected growth of the 2000 initial ice superparticles, (Rqj + Aa), 
as a function of initial radius Roj. Comparing the expected final 
particle radii with the actual final radii shows that for small par- 
ticle sizes the Monte Carlo method works fairly well, whereas 
for larger particle sizes most particles experience too little (zero) 
growth and a few grow several orders of magnitude more than 
expected. This means that one must be careful not to draw con- 
clusions from the growth of individual particles. However, sta- 
tistically speaking, i.e. averagin g over many particles, the result 
can be considered to be correct (Zsom & Dullem ondl 120081) . 

5. Results 

The results are divided into two parts: runs including only the 
atmospheric ice line (Table [TJ runs 1 and 2) and runs including 
both the radial and atmospheric ice lines (Table[TJ runs 3 and 4). 
We assess how the growth depends on turbulence strength in run 
2 and 4, its dependence on atmospheric ice line position in runs 
1, 2 and 3 and the dependence on box size in run 4. In all runs 
particle collisions are ignored, an assumption that we discuss in 
Sectionl6H 

5.1. Atmospheric ice line 

We start by exploring the atmospheric ice line, not including the 
radial ice line. This means that the simulation domain is located 
just outside of the radial ice line, with a c older midplane and 
hotter outer layers (iDullemond & Dominiki |2004|) . 

The thermal atmosph eric ice line in a typical pro toplanetary 
disc is located at z > 3 //(Ch iang & Goldreichlll997l) . However, 
the decrease of pressure with height above the midplane leads to 
a narrow radial zone just outside of the radial ice line where the 



atmospheric ice line is located at z < 3 H. The pressure ice line 
in a vertically isothermal disc is shown in Fig. |4] 

The height of the atmospheric ice line above the midplane 
Zice is varied from z; ce = 0.4// to z; ce = 3 H. The lower limit is 
the lowest value for which not all ice particles immediately sub- 
limate. A lower value gives a distance between the midplane and 
the ice line which is comparable to the length a particle moves 
during one time step (see Eq. [T7t . making it possible for all 
ice particles to escape from the condensation region before any 
growth has taken place. Periodic boundary conditions are used 
in both the radial and vertical direction. 

Fig. [5] shows the time evolution of a simulation with the at- 
mospheric ice line at z; ce = 1 H. The distribution of red vapour 
and blue ice particles is shown together with the size distribution 
of ice particles, for t = 0, t = 200 and t = 2000 Q~ l . Vapour 
diffuse into the grey condensation region, condensing onto al- 
ready existing ice particles. As small ice particles tend to diffuse 
out of the condensation region and sublimate, whereas larger 
ones stay in the midplane and experience continued growth, the 
result is a narrow size distribution with decimeter-sized ice peb- 
bles residing in a midplane layer. 

5.1.1. Varying the atmospheric ice line position 

Simulations were run with the atmospheric ice line placed at 
different heights z; ce above the midplane. Decreasing Zi ce cor- 
responds in a physical protoplanetary disc to placing the simu- 
lation box radially closer in towards the radial ice line, where 
a hypothetical z K e = would correspond to the position of the 
radial ice line, and thus assessing the narrow radial zone where 
the pressure ice line dominates. 

In Fig. [6] the particle growth for different z; ce is shown. The 
mean ice particle size (R) in units of Ri is shown as a function 
of time in units of Q x . The curves show, from top to bottom, 
z ice = 0.4// (black), z ice = 0.6 H (violet), z ice = 1.0// (blue), 
Z; ce = 1 .4 // (green), z; ce = 1.8 //(red) andzj ce = 3.0 H (yellow). 
As is clear from the figure, particles grow to larger sizes as the 
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Fig. 4. Pressure ice line for a vertically isothermal disc. The full black line denotes the ice line, whereas the coloured dotted lines 
show the location of 1 H, 2 H and 3 H. As particles are preferably found within a few scale heights from the midplane, growth due 
to crossing of the atmospheric ice line is only important within less than 1 AU of the radial ice line. 
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Fig. 5. State of simulation with only the atmospheric ice line included, for t — 0, 200 Sand 2000 Q~\ from left to right. Upper 
panels: Distribution of blue ice particles and red vapour particles, with the grey area representing the condensation zone. The 
size of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with 
size for visibility. Lower panels: Evolution of the size distribution of ice particles at corresponding times. is the number of ice 
superparticles, and the size is given as R/R\. At the water ice line, /fj * 1.3 m. 
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Fig. 6. Mean ice particle size as a function of time for different positions of the atmospheric ice line. The different lines 
show the growth for different distances from the atmospheric ice line to the midplane, Zj ce . From top to bottom Zi ce = 
0.4 H, 0.6 H, I. OH, 1 .4 H,\&H and 3.0 H. Ice particles grow faster, and to larger sizes, the closer to the midplane the atmospheric 
ice line is. Full lines denote simulations with 10 000 particles, dashed lines 1000 particles and dotted lines 100 000 particles. 



ice line is moved closer to the midplane. This is effectively due 
to the fact that for ice lines closer to the midplane, i.e. a nar- 
rower condensation region, vapour particles can easier penetrate 
to the larger ice particles residing in the midplane. Fig. |2]shows 
how the particle scale height H v decreases with increasing par- 
ticle radius R. Small particles have a scale height comparable to 
the gas scale height ff PiSma n ~ H, whereas for larger particles 
^p.iarge ^ H. Placing the ice line at a height larger than 1 H is 
thus equivalent to saying that even the smallest particles tend to 
stay within the condensation zone. There is therefore very little 
material available for particle growth. The few particles that do 
move across the ice line and sublimate will mostly redistribute 
material among the particles at the border of the condensation 
zone. Growth in runs with z; ce > 1 H is both slow, and to mod- 
est particle sizes. For ice line positions of z; ce > 3 H the particle 
growth is on the order of AR ~ 10~ 4 i?i or less in 1000 Q -1 . For 
the larger extent of the disc, where z; ce > 3H, the atmospheric 
ice line is thus of little importance for particle growth by con- 
densation. 

In Fig.|6]we show runs with 1000, 10 000 and 100 000 parti- 
cles for each ice line position. The negligible difference between 
runs with 10 000 and 100 000 particles in combination with a 
much lower computational cost motivates the use of 10 000 par- 
ticles as a reference resolution. 

5.1 .2. Varying the turbulent a-value 

Simulations were run with different values of a in order to test 
the influence of turbulence strength on the results. In Fig. |7]the 
mean particle size (R) is shown as a function of the height of the 
ice line above the midplane zi ce , starting from an initial particle 
size of 10~ 4 Ri . The system is shown when it has approximately 
reached an equilibrium, at t — 5000 Q l for a = 10~ 2 , shown in 
blue, at t — 25 000 Q~ l for a = 10~ 3 , shown in red, and at t — 
1000006T 1 for a = 10 4 , shown in yellow. Full lines represent 
modelled values and dotted lines are analytical estimates. As can 
be seen, lowering the strength of turbulence gives less growth. 



This is to be expected as the main effect of decreasing the level of 
turbulence is to lower the effective particle scale height, set by an 
equilibrium between sedimentation and turbulent diffusion. As a 
comparison, dotted lines show the analytical relation between 
mean particle size and particle scale height in a system where 
vertical gravity and turbulent diffusion is dominating. The scale 
height of particles in terms of the gas scale height is 

Setting Hp ~ (l/c)zi ce gives an analytical estimate of the maxi- 
mum particle size achievable by condensation, shown in Fig. [7] 

as 

Or f = a{ ,/ m l2 -l|- (54) 

We chose c = 3, but even in this case the analytical expression 
underestimates the resulting R slightly. Nevertheless our simu- 
lations show that particles grow approximately to a size where 
their scale-height is 1/3 of the ice line height. 

5.2. Atmospheric and radial ice line 

Including both the atmospheric and radial ice lines is done by 
letting the simulation domain represent a region around the ra- 
dial ice line. This means that ice particles can sublimate both 
by moving away from the midplane, and by moving closer to 
the central star. For a realistic atmospheric ice line position at 
Zice ^ 3 H, sublimation by moving across the radial ice line is 
however much more important than by moving across the atmo- 
spheric ice line. The simulation domain is 12 H in the vertical di- 
rection. In the radial direction the box size is varied from 1 .5 H 
to 5.25 H, with the number of particles and grid cells changed 
accordingly in order to keep the particle density and the effec- 
tive resolution fixed. Periodic boundary conditions are used in 
the vertical direction, and reflective boundary conditions in the 
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Fig. 7. Mean ice particle size as a function of the position of the atmospheric ice line and the strength of the turbulence. The mean 
particle size (R) is shown as a function of the distance from the atmospheric ice line to the midplane z; ce , in gas scale heights H. The 
system is shown when the particle sizes have approximately reached an equilibrium, at t — 5000 Q l for a = 10~ 2 , shown in blue, at 
t = 25 000 Q l for a = 10~ 3 , shown in red, and at t = 100 000 Q l for a = 10~ 4 , shown in yellow. The initial particle size is 10 4 i?,, 
corresponding to approximately 0.1mm at the ice line. Solid lines are modelled values and dotted lines are analytical values, set by 
Eq.[53] The modelled and analytical values follow the same slope in the applicable range. 



1.000 




1000 



Fig. 8. Mean ice particle size (R) as a function of time in in- 
verse Keplerian orbits Q~ l for simulations including both radial 
and atmospheric ice lines. The height of the atmospheric ice line 
above the midplane Zke is denoted as: black 0.4 H, violet 0.6 H, 
blue 1.0 H, green 1.4 H, red 1.8 H and yellow 3.0 H. For the 
runs with Zi ce = 3.0 H, the dashed line corresponds to 1000 par- 
ticles and the dotted line to 100000 particles. Full lines denote 
runs with 10000 particles. The position of the radial ice line is 
Ice = _ 0.3 H, and the strength of turbulence is a = 10~ 2 . 



radial direction. In the vertical direction there is in effect no dif- 
ference between +z and — z, as the conditions above the midplane 
mirror those below it. In the radial direction, when including the 



radial ice line, there is no such symmetry. At r — r m i„ conditions 
such that ice sublimates prevail at all z, whereas at r = r max the 
midplane is part of the condensation zone, where ice particles 
can exist without sublimating. Using periodic boundary condi- 
tions in the radial direction would therefore artificially introduce 
vapour in the system, causing a too large particle growth. 

5.2.1. Varying the atmospheric ice line position 

Including the radial ice line greatly increases the growth ef- 
ficiency As shown in Sec. I5.ll models with only the atmo- 
spheric ice line gives negligible growth for ice line positions 
of Zice ^ 3 H. However, taking also the radial ice line into ac- 
count leads to particle growth beyond 0.1 R i, corresponding to 
decimeter-sized pebbles at the ice line, within 10006T 1 . 

Fig.[8]shows how the mean particle size (R) as a function of 
time varies with atmospheric ice line position z; ce when also the 
radial ice line is included. The colours show different heights of 
the atmospheric ice line above the midplane, with black denot- 
ing Zke = 0.4//, violet Zi ce = 0.6 H, blue Zj ce = 1.0//, green 
Zke = 1.4H, red Zi ce = 1.8 /f and yellow Zj ce = 3.0 H. Growth 
is somewhat faster for Zi ce close to the midplane and slower for 
Zice further away from it. As previously discussed, this is due 
to the larger supply of material available for growth in a nar- 
row condensation region, as compared to a wider condensation 
region. For ice line position above 3 H, a similar growth as for 
Zke = 3 H is found. 

Growth stops at (R) ss R\ for all Zke- This is due to radial 
drift towards the star, which eventually removes all ice particles 
from the simulation. From Eq. [30] it can be seen that the radial 
drift inwards peaks for particles of (R) « R\. As particles reach 
this size they thus drift inwards, past the radial ice line, and sub- 
limate. 

We also show the growth for ice line position of Zi ce = 3.0 H 
for a smaller (n p = 1000) and a larger (n p = 100000) number 
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Fig. 9. State of an extended simulation box with both the radial and the atmospheric ice line included, for t — 0, 100 and 1000 Q -1 , 
from top to bottom. The grey area represents the condensation zone, and ice and vapour is shown in blue and red, respectively. The 
sizes of the blue dots are proportional to the size of the ice particles and the number of particles shown is inversely scaled with size 
for visibility. The turbulent <?-value is 10~ 2 . 



of particles in Fig. [8] in addition to n v = 10000 particles. For 
particle sizes R < 0.1 Ri we find converging results, whereas 
for larger particle sizes the growth is dominated by statistical 
fluctuations caused by runaway-growth of a few lucky particles. 
Our condensation model is constructed to give correct results 
in a regime where many particles compete for vapour, i.e. up to 
pebbles of a few decimeters, but may be unphysical in the regime 
of meter-sized boulders and beyond. 

5.2.2. Varying the radial extent of the simulation and the 
turbulent a-value 

To assess the importance of the size of the simulation domain, 
we ran simulations where the box size was varied in the radial 
direction. Fig. [9] shows the time evolution of a run where the 
radial extent of the simulation box has been increased, so that 
fmax = 5.25 H. The most significant growth takes place close to 
the radial ice line, but due to radial mixing large particles can be 
found throughout the entire radial extent of the simulation box. 
In Fig. [10] the partial growth in four equally sized subdomains 
can be compared to the total growth in the whole simulation do- 
main, as shown in Fig. [9] This illustrates that the total particle 
growth is dominated by the growth near to the radial ice line. 

The left panel of Fig. II II shows the particle growth for dif- 
ferent radial extents of simulation domains for a strongly turbu- 
lent disc, a = 10~ 2 . The growth timescale increases with larger 
simulation domains, as the amount of ice particles is increased 
while the supply of water vapour across the ice line remains un- 
changed. Results from simulations with weak turbulence are pre- 
sented in the middle (a = 10~ 3 ) and right (a = 10~ 4 ) panels of 
Fig. QT| These low-turbulence simulations show that, although 
slowly, particles can grow via condensation also in dead zones. 

6. Discussion 

Our results show that ice condensation is an efficient way to form 
centimeter-sized and decimeter-sized pebbles. Ice condensation 
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Fig. 10. Growth in different radial subdomains of the simulation 
box shown in Fig. [9] Black shows the average particle size in the 
total box, whereas yellow, red, green and blue gives the average 
particle sizes in the different subdomains, from left to right. The 
total particle size at any time is dominated by the growth near 
the radial ice line. 



does not suffer from bouncing and fragmentation barriers and 
could be the dominant mode of growth in protoplanetary discs, 
to sizes where particles can concentrate strongly in the turbulent 
gas and continue towards planetesimal sizes by gravitational col- 
lapse. However, our results rely on a number of assumptions and 
simplifications which we discuss here. 
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Fig. 11. Particle growth for different radial extents of the simulation box and turbulence strengths. The mean particle size (R) is 
shown as a function of time in inverse Keplerian orbits Q~ l for a = 10~ 2 (left panel), a = 10~ 3 (middle panel) and a - 10~ 4 (right 
panel). The size of the box has been extended in the radial direction, with r max colour coded as: black 0.75 H, blue 2.25 H, red 3.75 
H and yellow 5.25 H. For all runs Zmm = -6.0 H, z max = 6.0 H and r m [ n = -0.75 H. The inlay figure shows the simulation domain, 
with the considered radial extents marked. 



6. 1 . Other ice lines 

This model has been developed primarily to investigate the wa- 
ter ice line at r ~ 3 AU. Similar ice lines, or condensation fronts, 
exist also for other chemical species, and as this model is scal- 
able to any r it can easily be adapted to include any condensation 
front. Of importance for planet formation purposes are in partic- 
ular the ammonia, methane, carbon monoxide and molecular ni- 
trogen ice lines further out in the protoplanetary disc, and the nu- 
merous silicate condensation fronts further in towards the central 
star. A global disc model, including all condensation fronts of the 
most important chemical species in the protoplanetary disc, will 
be an important extension of this work. 



6.2. Disc evolution 

The systematic motion of gas accretion onto the central star has 
been ignored in this work. However, compared to the time scales 
of particle growth by condens ation the accretion proc ess is slow 
enough to be safely neglected (Hartm ann et all fl 998). 

The ice line has throughout this work been considered as 
fixed. Seen over the life time of the disc this is not true, as 
the ice line probably moves both in a systemati c way over lon g 
time scales and due to shorter heating events (Armitagel |201 ll) . 
However, as the time scale for growth by condensation found in 
this work is on the order of r c w 1000 Q , which is very short 
compared to the life time of the disc Tdi sc ~ 3 Myr, the assump- 
tion of a fixed ice line is reasonable. 

We also ignore the potential effect the processes of con- 
densation and sublimation have on the ice line position, 
since the resulting release and absorption of energy is rela- 
tively small compared to th e total thermal energy of the disc 
dStevenson & Lunindfl988h . 



6.3. Ice nuclei 

We have modelled ice particles as one-component particles; ho- 
mogeneous water ice particles. In reality these particles would 
have (at least) a two-component composition, that of a rocky 
core and an ice mantle. This is due to the fact that supersaturated 
water vapour is needed for homogeneous nucleation, i.e. for the 
spontaneous formation of a new ice particle. As this is typically 
not the case in the conditions prevailing in a protoplanetary disc, 
water vapour instead condenses onto an existing grain, either a 
rocky dust particle, or a particle with an already existing ice layer 
on top of the dust grain. Whether vapour preferably condenses 
onto a bare dust grain or an ice particle depends on material 
properties and on their respective sizes. 

Ignoring the small refractory dust particles present in the disc 
possibly leads to an over-estimation of the growth of the larger 
ice particles. For a fixed amount of mass, smaller dust particles 
have a larger combined surface area than larger particles, and 
therefore there is a higher probability that a vapour particle con- 
denses onto a smaller dust superparticle than onto a larger su- 
perparticle. This is for the case when the material properties of 
dust and ice are such that two equally large dust and ice particles 
have equal probabilities of having a vapour particle condensing 
onto them. Instead of growth of a few large ice particles that sed- 
iment towards the midplane and thereby are protected against 
sublimation, there would thus be a continuous sublimation and 
condensation of the small ice particles at the ice line. This ef- 
fect can already be seen in our condensation model, as small 
particles are more likely to cross the ice line and sublimate, and 
conversely, to grow by condensation, but could be amplified by 
the introduction of refractory grains. 

Material effects might make this problem more severe. The 
dust grain composition is assumed to be similar to that of inter- 
stellar dust grains, either being pure interstellar dust grains or 
conglomerates thereof. This gives us the most important dust 
species as silicates and carbonaceous material , with s ilicate s 
typically assumed to be the dominating one dPrainel 120031) . 
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Silicates have material properties that makes t hem more effi- 
cient as ice nuclei than pure water ice particles dGoumans et all 
2009), implying that the effect of small particle growth at the 
expense of the larger ones is amplified when taking the material 
properties into account. 

There are however a number of possible solutions to this 
problem, in terms of reasons for the vapour to condense onto ice 
particles instead of the small dust grains. Firstly, the dust grains 
are not likely to all have the same size. Rather, they would have 
a size distribution w here the smallest grains can be nanometer- 
sized dDraineir2003l) . The largest dust particles of relevance are 
the ones of the largest size that are strongly enough coupled 
to the gas to be present at the relevant vertical distance from 
the midplane. This is around micrometer-scale (see Fig. |3}. The 
Kelvin curvature effect states that the equilibrium vapour density 
for a curved surface is higher than for a flat surface, meaning that 
vapour more easily condenses onto a fiat surface than a curved 
one, an effect which is m ost important f or the smallest particles, 
around nanometer-sizes dSironoi 1201 ]]). Therefore vapour will 
effectively prefer condensing onto the larger particles, leaving 
the smallest ones as bare dust grains. 

Secondly, small grains are likely to be hotter than larger par- 
ticles, moving the ice line of small particles closer to the mid- 
plane than the ice line of large particles. The temperature of the 
grains is set by the balance between absorption and emission ef- 
ficiency. As grains absorb and emit radiation most efficiently in 
wavelengths smal ler than and up to their own size, there is a size- 
dependent effect dKrivov et all l2008t IVitense et all 120121) . The 
spectral energy distribution of a solar-type star peaks at wave- 
lengths of about 500 nm. All grains of this size and larger there- 
fore absorbs incoming radiation with approximately equal ef- 
ficiency. However, grains emit in infra-red wavelengths (larger 
than 1 //m), and thus the larger grains can emit more energy 
than the smaller ones. The resulting temperature for the smaller 
grains is therefore higher than for the larger grains. Also the ma- 
terial affects the temperature of the grains. As water ice is nearly 
transparent in visible light, where a solar-type star peaks, an ice 
particle is heated less by t he stellar radiation than a particle of a 
more absorbing material dLecar et all H0Q6). Both mechanisms 
are valid for grains residing in the atmosphere of the disc, where 
the disc is optically thin to the stellar radiation, but not for grains 
in the midplane, where the received stellar radiation is mainly the 
re-emitted infrared radiation from neighbouring grains. 

Further, there is a possibility that material effects might pre- 
vent water vapour from condensing onto the dust grains. As 
mentioned above, bare silicate grains are more efficient as ice 
nuclei than ice-cove red ones. Howev er, for carbonaceous grains 
the opposite is true (Papoulaij |2005l) . meaning that vapour con- 
denses more easily onto a water ice particle than onto a bare car- 
bonaceous grain. For a grain population with ice coated particles 
and bare carbonaceous grains, the carbonaceous grain would 
therefore be left bare, whereas the ice grains would continue 
growing. 

6.4. Particle collisions 

We have neglected collisions between particles throughout this 
work. Depending on material properties and relative velocities of 
the particles, collisions can lead to increased growth via coagu- 
lation, to fragmentation or to bouncing. In general, small parti- 
cles tend to stick together to a larger extent than large particles, if 
colliding, whereas larger particles are more prone to bouncing or 
fragmenting. Collisions amongst large particles therefore have a 
tendency to be destructive, whereas small-particle collisions are 



more likely to favour growth (or at least not counteract it). For 
silicate particles growth via collisions involving equal-sized par- 
ticles is very difficult beyond millimeters, however collisional 
growth involving small parti cles colliding with a large target 
is possible, although slow dWurm et all 120051; iJohansen et all 
l2008l:IWindmark et alll2012h . 

Whereas extensive experimental data exists on collisions be- 
tween rocky particles, th e outcome of i ce p article collisions is 
not equally well-known. iBridges et all (1 19961) performed low- 
velocity collisions (v le i < 5cms _1 ) between ice pebbles and 
found that ice particles covered with a frost layer have an in- 
creased stickiness compared to rocky particles. A mechanism of 
collisional fusion, in which ice particles undergo a phase change 
during collision, has been suggested as a way for particles to 
stick also i n collisions with higher velocities (1ms 1 < v re i < 
100ms 1 ) dWettlauferll20Tol) . 

Collision experiment for higher velocities have been used 
to derive criticia l relative v e locitie s above which fragmentation 
can occur, v CI i t . Iffiga et all d 19961) found a critical velocity of 
v re i » lms' 1 , but other groups have found significantly larger 
values (Arakawa, 1999; lArakawa et all 120021) . Computer sim- 
ulations suggest critical velocities of up to v re i ~ 100ms 1 
dBenz& Asphaugill999l) . 

Collision velocities for millimeter-sized particles are set by 
the turbulent velocities and are expected to be on the order of a 
few m s _1 , whereas collisions involving larg er particles approach 
the radial drift velocity, v re i m 50ms _1 dBrauer et al. , 2008). 
Both coagulation and fragmentation are therefore expected to 
occur as a result of collisions. A likely scenario is that collisions 
between large particles lead to fragmentation, where the result- 
ing small ice fragments are later swept up in collisions with other 
particles. 

An important implication of collisions is that it provides a 
natural means of removing small dust grains released from the 
ice particles when sublimating. As these dust grains are very effi- 
cient ice nuclei, they might prevent growth of already large par- 
ticles by "stealing" all the vapour, as discussed in Section 1631 
However, the small dust grains are likely to be swept up by larger 
ice particles in collisions. This does not increase the growth sig- 
nificantly, but has the important benefit that the small dust grains 
are removed so that vapour has to condense onto growing ice 
particles. 

We plan to include the effects of multiple condensable 
species, refractory ice nuclei and particle collisions in future 
work. The present work highlights that simple turbulent dy- 
namics can cause significant particle growth by condensation of 
volatiles, motivating follow-up studies with increasingly realis- 
tic models for the condensation and collision processes. 

7. Conclusions 

In this paper we have shown that ice condensation is an impor- 
tant particle growth mechanism that must to be taken into ac- 
count in models of early planet formation. As the more thor- 
oughly investigated mechanism of coagulation is inefficient in 
forming particles larger than centimeters, growth by condensa- 
tion is an important mechanism that could complement coagula- 
tion or even be the dominant particle growth mode. 

Our results show that growth by condensation near ice lines 
is rapid and results in large particle sizes. The growth time scale 
from millimeter-sized dust grains to decimeter-sized pebbles is 
t c w 1000 years. Significant growth is obtained for a range of 
turbulent a-values from 10~ 4 to 10~ 2 , where the higher value 
corresponds to the strength of turbulence that has been inferred 
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from observations of protoplanetary discs, and the lower value 
has been suggested in a dead zone with weak turbulence stirred 
by the narrow active layers. 

An implication of our model is a lower column density of 
vapour in the entire region outside the radial ice line compared 
to the inner part of the disc, caused by condensation onto ex- 
isting grains and subsequent sedimentat ion of particles . This is 
in agreement with observations of CO (lOi et al.L 1201 ll) . where 
a drop in abundance was found outside of the CO ice line. 
Similarly, the outer disc regions hav e been observationally in - 
ferred to be depleted in water vapour (Hogerheii de et al.L |2011|) . 
Observed remnant water and CO vapour in the outer disc shows 
the importance of the coexistence of a cold midplane and an up- 
per atmospheric layer where vapour can still exist. 

Ice condensation can also explain obser vations of large quan - 
tities of pebbles in protoplanetary discs (Wilner et al., 2005). 
Such pebbles are crucial in planet formation models, as they are 
the preferred starting size for planetesimal formation by clump- 
ing via streaming instabilities, in pressure bumps and in vortices, 
followed by gravitational collapse. Once large planetesimals 
have formed, continued pebble accretion is a v ery efficient path 
to formation of the cores of gas and ice giants (lOrmel & Klahr, 
20 lOt lLambrechts & Johanseni 1201 2t iMorbidelli & Nesvornvl 
201% . 

Ice condensation is a natural consequence oftemperature gra- 
dients in protoplanetary discs. With our simulations we have 
shown that condensation is an efficient growth mechanism with 
the potential to explain the formation of decimeter-sized peb- 
bles in protoplanetary discs, thereby providing the missing link 
to further growth into planetesimals and planets. 
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